# 主函数

# 预处理
data_process <- function(data){
  # list to frame
  df = as.data.frame(data)
  # 获取指定列
  df_new = df[, c(1, 3)]
  
  # 列名重置
  colnames(df_new)[2] <- "freight_volume"
  
  # 日期格式转换
  df_new$month <- as.Date(paste0(df_new$month, 1), "%Y年%m月%d")
  start_date = min(df_new$month)
  end_date = max(df_new$month)
  # 排序
  df_new <- df_new[order(df_new$month),]
  
  # frame to ts 格式转换
  series <- ts(df_new[[2]], 
               start = c(format(start_date, "%Y"), format(start_date, "%m")), 
               end = c(format(end_date, "%Y"), format(end_date, "%m")),
               frequency = 12)
  # print(series)
  
  ###  切片，分训练集和测试集
  series_train = window(series, start = 2015, end = 2019+6/12)
  series_test = window(series, start = 2019+8/12)
  # print(series_train)
  # print(series_test)
  return(list(series_train, series_test))
}


# 代码需要修改的路径
# common_path <- "D:/Learn/Allen/code/gitee_project/times_series"
common_path <- "F:/code/time_series"
work_path <- paste(common_path, "code/R", sep="/")
file <- paste(common_path, "data/运输量.xlsx", sep="/")

# 设置工作路径
setwd(work_path)

# load xlsx
source("read_file.R")

sheet_data <- read_xlsx_data(file)


# 预处理
data_info = data_process(data = sheet_data)
series_train = data_info[[1]]
series_test = data_info[[2]]

# 模型 & 参数选择
# ARIMA
fit_arima <- auto.arima(series_train)
checkresiduals(fit_arima)

# ETS
fit_ets <- ets(series_train)
checkresiduals(fit_ets)

# 预测效果
evalute_arima <- fit_arima %>% 
  forecast(h = length(series_test)) %>%
  accuracy(series_test)

evalute_arima[, c("RMSE", "MAE", "MAPE", "MASE")]

evalute_ets <- fit_ets %>% 
  forecast(h = length(series_test)) %>%
  accuracy(series_test)

evalute_ets[, c("RMSE", "MAE", "MAPE", "MASE")]

######### 动态谐波回归：傅里叶变换
library(fpp2)

plots <- c()
for (i in seq(6)){
  xreg_train <- fourier(series_train, K=i)
  xreg_test <- fourier(series_test, K=i, h=length(series_test))
  
  fit <- auto.arima(series_train, 
                    xreg = xreg_train, 
                    seasonal = FALSE, 
                    lambda = 0)
  
  plots[[i]] <- autoplot(forecast(fit, xreg=xreg_test)) + 
    xlab(paste("K=", i, " AICC=", round(fit$aicc, 2)))
  
}

library(gridExtra)

grid.arrange(grobs=plots, ncol = 2)

